Effective Temperatures of a Driven System Near Jamming 
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Fluctuations in a model of a sheared, zero-temperature foam are studied numerically. Five differ- 
ent quantities that reduce to the true temperature in an equilibrium thermal system are calculated. 
All five have the same shear-rate dependence, and three have the same value. Near the onset of 
jamming, the relaxation time is the same function of these three temperatures in the sheared system 
as of the true temperature in an unsheared system. These results imply that statistical mechanics 
is useful for the system and provide strong support for the concept of jamming. 
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Statistical mechanics describes the connection between 
microscopic properties and collective many-body proper- 
ties in systems in thermal equilibrium. There is no equiv- 
alent formalism for driven, athermal systems. Nonethe- 
less, recent phenomenological approaches [Ql that assume 
thermal behavior are surprisingly successful in describ- 
ing driven glassy materials such as sheared foam. Foam 
is a dense packing of bubbles in a small amount of liq- 
uid, and is athermal because the thermal energy is much 
smaller than the typical energy barrier for bubbles to 
change their relative positions As a result, quies- 
cent foam is jammed it is disordered and has a yield 
stress. If foam is steadily sheared, however, it is pushed 
over energy barriers and flows as different bubble pack- 
ings are explored. However, it is unclear if this degree of 
ergodicity is enough to lead to thermal behavior. 

In this Letter, we test the assumption that a sheared 
foam can be modeled as a thermal system with a tem- 
perature that depends on shear rate. We conduct numer- 
ical simulations of a simple model of sheared foam and 
measure five quantities that all reduce to the true tem- 
perature in a thermal system. Although these quantities 
must all have the same value in an equilibrium thermal 
system, there is no guarantee that they should be the 
same in the steadily-sheared model foam. Remarkably, 
three of the effective temperatures are the same and all 
five have the same shear-rate dependence. Our results for 
four of the effective temperatures are shown as a function 
of shear rate in Fig. ^for two different size distributions of 
the bubbles. In a companion paper, Berthier and Barrat 
reach similar conclusions for a sheared thermal system 
These results suggest that statistical mechanics is 
indeed useful for describing driven jamming systems. 

Our bubble dynamics (BD) simulations are carried out 
in two dimensions on Durian's model of foam Q. The 
bubbles are circles with diameters assigned from one of 
two different diameter distributions. The first (polydis- 
perse) distribution is flat from 0.2 to 1.8 times the av- 
erage bubble diameter and is zero otherwise. The sec- 



ond (bidisperse) distribution consists of equal numbers 
of small and large bubbles of diameter ratio 1.4. Pairs of 
bubbles only interact, via a repulsive spring, when they 
overlap; this approximates the energy cost of bubble de- 
formation [Q . There is also a frictional force proportional 
to the velocity difference between a bubble and the av- 
erage flow at its position. The system is fully periodic, 
with flow in the x direction and a shear gradient in the 
ij direction imposed using the Lees-Edwards boundary 
condition js]. 
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FIG. 1. Four effective temperatures, Tp (circles), Txy (up- 
ward triangles), Te (squares) and T/ (downward triangles) 
calculated as a function of shear rate. Tj has been rescaled 
by 4 and 2.7 for the polydisperse and bidisperse systems, re- 
spectively, to collapse on the others. Insets: Tp vs. 7 for 
different system sizes. 
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We scale lengths by the average diameter d, energies 
by kd^ where k is the spring constant, and time scales 
by To = b/k, where b is the friction coefficient. Thus, 
the dimensionless shear rate 7 is the Deborah number. 
Unless otherwise specified, the systems contain N — 400 
bubbles at an area fraction of = 0.9 (well above ran- 
dom close-packing at 0.84). Averages are typically taken 
over time snapshots separated by 0.1 in strain over a to- 
tal strain of 10 for 100 (polydisperse) or 2 (bidisperse) 
different initial configurations. 
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FIG. 2. (a) Shear stress autocorrelation functions at differ- 
ent shear rates for the polydisperse system (see Def. 2). (b) 
Variance of energy fluctuations plotted against average en- 
ergy squared (Def 3). The lines are fits to {{SEf) = m{Ef, 
where m = 3.9 x 10~^ (polydisperse). (c) Force distributions 
P{f) at different shear rates for the polydisperse system. The 
tail of P{f) is fitted to a Gaussian centered at zero to obtain 
the effective temperature Tf (Def. 4). 

Three of the 5 calculated effective temperatures 
are based on linear response relations or fluctuation- 
dissipation relations. Similar definitions have been used 
to characterize non-equilibrium systems in the contexts 
of weak turbulence Q , aging of glassy systems gran- 
ular packings ||,^ and sheared aging systems ; such 
effective temperatures can control the direction of heat 



flow and thus play the role of temperature in the thermo- 
dynamical sense in certain nonequilibrium systems with 
small energy flows pT| ]. 

Definition 1: pressure fluctuations. In an equilibrium 
system at fixed N, T and area A, the variance of the 
pressure is given by [|| 



(1) 



where = -A{d{p)/dA)T is the inverse isothermal 
compressibility, A is the area of the system, p is the pres- 
sure, X is the hyper- virial |^ , and the Boltzmann constant 
is unity. In our driven, athermal system we thus define 



A((V)) 



(2) 



so that Tp reduces to the true temperature in an equi- 
librium thermal system. To measure the compressibility, 
we perturb the system area A and measure the resulting 
value of (p) during shear to calculate the derivative. In 
simulations of a quiescent system in thermal equilibrium 
with the same potential and polydispersity, we find that 
with comparable statistics, Eq. ^ yields results within 5% 
of the simulation temperature. 

Definition 2: shear stress fluctuations. The viscosity 
of an equilibrium system is related to the integral over 
the shear stress autocorrelation function Q : 



A f°° 

= - dt{a^y{t)a^y{0))c 



(3) 



Again, we use Eq. |3|to define T^y, calculating the steady- 
state shear viscosity from rj — {(Jxy)/'y. The stress auto- 
correlation function is shown for several different shear 
rates in Fig. ||(a). With decreasing 7, the correlation 
time increases approximately linearly, while the variance 
{{ScTxy)'^) decreases and then saturates. 

Definition 3: energy fluctuations. The constant- 
volume heat capacity of an equilibrium system is related 
to energy fluctuations: 



d{E) USE)') 



dT 



rp2 



(4) 



This can be rearranged and integrated on both sides to 
provide a definition of T To calculate Te we must 
extract the relation between the variance of the energy 
fluctuations, ((5£')^), and the average energy, (E). The 
results are shown in Fig. |^(b) . We find that the variance 
scales as the square of the average energy, as shown by 
the line-fits to the data. The fits imply that (E) oc Te, 
as can be seen by substitution into Eq. |^. In fact, we 
find {E) — O.GANTe for the polydisperse system and 
(E) — 0.56NTe for the bidisperse system, where N is 
the number of bubbles in our two-dimensional system. 
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These results resemble equipartition, except that the co- 
efficient of NT is not unity. Our potential is a harmonic 
repulsion with finite range, so equipartition is not exact. 

Definition 4- force distribution. The temperature can 
be extracted from the tail of the distribution of inter- 
particle normal forces P{f) |13|. The force distribu- 



tion is directly related to the pair correlation function 
g{r) = ex.p{—f3u{r))y{r) in a system with a pair poten- 
tial u{r). For r sufficiently small, the exponential term is 
a much stronger function of r (and hence of /) than y{r). 
As a result, the tail of P{f) depends on the interparticle 
potential and the temperature: 



P(/) « exp(-/V2r). 



(5) 



We fit the tail to extract Tf. Eq. || applies only to a 
monodisperse system. For a bidisperse system, we mea- 
sure the distribution for small particles interacting with 
small particles only. For the polydisperse case, however, 
we use the force distribution for all particles (Fig. |^(c)). 
From MD simulations of an equilibrium system with the 
same potential, we have verified that this assumption 
should lead to no more than a 30% error in the com- 
puted Tf. 
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FIG. 3. The Stokes-Einstein relation. The circles represent 
the diffusion constant D, measured by integrating the veloc- 
ity autocorrelation function. The triangles represent Tp/Crjd, 
where C is a constant chosen to obtain the best fit with D. 

Definition 5: Stokes-Einstein relation. In equilibrium, 
the diffusion constant D satisfies 



D 



T 
Chid' 



(6) 



where C depends on system size in two dimensions. Thus, 
Eq. |6| defines a temperature To up to an unknown con- 
stant. Here we use fixed boundary conditions in the y- 
direction We measure D in the y-direction in two 

different ways and find good agreement: we integrate the 
velocity autocorrelation function, and we measure the 
displacement distribution as a function of time from an 
initial starting position. It is difficult to extract To be- 
cause D and rj vary by several orders of magnitude over 



the range of 7 studied, while their product varies by less 
than an order of magnitude and has a lot of scatter and 
large error bars. This definition is therefore most use- 
ful as a consistency check: we use Tp from Def. 1 and 
vary C to obtain the best agreement between the left and 
right sides of Eq. |^. Fig. ^ shows that the Stokes-Einstein 
relation is indeed obeyed. 

Our results for the other four effective temperatures 
(Defs. 1-4) are plotted in Fig. |l|. The insets show that 
the temperatures do not depend on system size for N suf- 
ficiently large. All the temperatures have the same shear 
rate dependence over 4 decades of 7. Since we do not 
have a first-principles calculation of C, the magnitude 
of Td is not known. However, Tf is different in magni- 
tude (but not in 7-dependence) from the remaining three. 
One possible reason for the discrepancy is that Tf mea- 
sures the properties of the tail of P{f), while Tp, T^y and 
Te measure fluctuations around average quantities. This 
suggests that the underlying probability distribution may 
not be a Boltzmann distribution. 

Fig. suggests that the effective temperatures ap- 
proach nonzero constants T^^g in the limit 7 — > 0. In 
fact, this is expected from their definitions. As long as 
bubbles overlap as 7 — > 0, the force distribution P{f) is 
nonzero for / > 0, and yields T^ > 0. The other tem- 
peratures should also be nonzero; in the zero shear rate 
limit, the angular brackets indicate configurational aver- 
ages rather than time averages. This suggests an inter- 
pretation of the limiting value; T^f^ should correspond to 
the glass transition temperature Tg. By shearing the sys- 
tem and calculating averages over times long compared 
to the relaxation time (which scales as I/7), we are de- 
manding that the system is ergodic. Therefore, the only 
temperatures accessible to us are above Tg. 

We have checked the interpretation of T^ff as Tg by 
conducting equilibrium molecular dynamics (MD) simu- 
lations on a system with the same interaction potential, 
packing fraction, and size distribution (bidisperse) as in 
our bubble dynamics (BD) runs. It has been proposed 
that jamming systems such as this one can be described 
by a phase diagram sketched in the inset to Fig. IJ. 
This diagram shows that jamming occurs as T is low- 
ered, the packing fraction (j) is raised, or the applied shear 
stress axy is lowered, and has been shown to be a useful 
way to represent experimental data jl^ . The BD results 
in Fig. ^ correspond to the trajectory marked "BD" in 
the inset to Fig. ^. In the MD simulations, we have re- 
moved the frictional term from the equations of motion 
and added the incrtial term and true temperature so as 
to approach jamming along the trajectory marked "MD." 
In Fig. ^, we show the results for the relaxation time t 
as a function of T from MD and BD, at two different 
packing fractions, 4> = 0.85 and (j> = 0.90. The lower 
packing fraction is just above random close-packing. In 
MD, we measure r from the decay of the intermediate 
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scattering function |T^]. In BD, we plot t = c/7, where 
c is a constant chosen to best fit the MD data [0 , as a 
function of Te- 

Fig. ^ shows that the dynamics are the same for a sys- 
tem approaching jamming by two different trajectories: 
decreasing T and decreasing 7 (or equivalently, decreas- 
ing {oxy))- Thus, Tg from MD at zero apphed shear is 
the same as T^ff from BD. A similar result was found 
previously for a sheared thermal Lennard-Jones mixture 
[0. Fig. H also shows that the functional form of the 
slowing down of the dynamics is the same along both 
trajectories for low T. This is consistent with previous 
work showing that the dynamics of the sheared model 
foam can be described by a Vogel-Fulcher form An 
attempt to collapse the data for the two different us- 
ing an Angell fragility plot shows that the fragility 
depends on packing fraction. 



We speculate that the concept of effective temperature 
is useful if nearly all the power supplied by the driving 
force is dissipated by fluctuations. In the model studied 
here, all of the power is, by construction, dissipated by 
fluctuations-the frictional force is proportional to the dif- 
ference between the velocity of a bubble and the average 
shear. In systems in the stick-slip regime near jamming, 
fluctuations typically are large compared to the average 
flow This suggests that the concept of effective tem- 
perature should be useful for any system near the onset 
of jamming. 
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NSF-DMR-0087349 (IKO,CSOH,AJL) and NSF-DMR- 
9722646 (CSOH,SRN). 
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FIG. 4. An Arrhenius plot of the relaxation time as a 
function of temperature from MD (gray symbols) and BD 
(black symbols) simulations at two different packing frac- 
tions, 4> ~ 0.85 and (p — 0.90. Inset: the jamming phase 
diagram. Our bubble dynamics simulations follow the trajec- 
tory marked BD, while our molecular dynamics simulations 
follow the trajectory marked MD. 

The results in Fig. |^ provide confirmation of the idea 
underlying the jamming phase diagram, namely, that the 
jammed region controls the behavior nearby so the dy- 
namics should not depend on the direction along which 
the jammed region is approached p^,p[. It also suggests 
that one can collapse the cr— axis onto the T— axis us- 
ing Toff, and that a system will jam once fluctuations, 
whether thermal or shear-induced, are sufficiently small. 

Our main finding that 5 different effective tempera- 
tures have the same 7 dependence raises the need for a 
criterion for when the concept of effective temperature 
might be useful. We suggest such a criterion based on 
the idea underlying the fluctuation-dissipation relation. 
An analogous concept applies to a steady-state driven 
system, because the average power supplied to the sys- 
tem must be balanced by the average power dissipated. 
The power can be dissipated in two ways-by the aver- 
age flow and by fluctuations around the average flow. 
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